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ABSTRACT 

We derive the evolution equations describing a thin axisymmetric disk of gas and stars with an arbi- 
trary rotation curve that is kept in a state of marginal gravitational instability and energy equilibrium 
due to the balance between energy released by accretion and energy lost due to decay of turbulence. 
Rather than adopt a parameterized a prescription, we instead use the condition of marginal gravita- 
tional instability to self-consistently determine the position- and time-dependent transport rates. We 
show that there is a steady-state configuration for disks dominated by gravitational instability, and 
that this steady state persists even when star formation is taken into account if the accretion rate 
is sufficiently large. For disks in this state we analytically determine the velocity dispersion, surface 
density, and rates of mass and angular momentum transport as a function of the gas mass fraction, 
the rotation curve, and the rate of external accretion onto the disk edge. We show that disks that are 
initially out of steady state will evolve into it on the viscous timescale of the disk, which is comparable 
to the orbital period if the accretion rate is high. Finally, we discuss the implications of these results 
for the structure of disks in a broad range of environments, including high redshift galaxies, the outer 
gaseous disks of local galaxies, and accretion disks around protostars. 

Subject headings: accretion, accretion disks — galaxies: evolution — galaxies: ISM — instabilities — 
ISM: kinematics and dynamics — turbulence 
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1. INTRODUCTION 

In the past few years, observational and theoretical ad- 
vances in many areas have led to intense study of the role 
of accretion and gravitational instability in determining 
the structure and rates of transport through disks. In 
the high redshift universe, clum py star-forming galax- 



ies a t redshifts z ~ 2 — 3 (e.g. lElmegreen et al. 
20051: iGenzel et all 120081: iForster Schreiber et alJ 



2004, 



2009; 



Cresci et al.ll2009l) appear to be undergoing rapid accre 



tion, and also have velocity dispersions that are much 
larger than those p resent in local galaxies. Numeri 
cal si mulations (e.g. iDekel et al.ll2009al iCeverino et al 
20091: lAgertz et all 20091 : iBournaud fc Elmegreenl 1200 



McNallv et al.l 12009 ) suggest that the large velocity dis 



persion and the massive clump morphology are both 
produced by a combination of gravitational instability 
and rapid external accretion. Around active galactic nu- 
clei, radiative cooling pushes thin accretion disks into a 
state of gravita tional instability (jShlosman 
IGoodmanl 120031 ) . and in this state their accretion rates 
and structures ma y be determin ed by gravitationally- 
driven turbulence ammie 2001). Closer to home, ac- 
cretion disks around protostars of mass > 1 Mq are ex- 
pected to experience strong gr avitational instabilit y for 
a significant part of their lives (jKratter et al.|[2008h due 
to a combination of rapid accretion and strong radiative 
cooling. Numerical simulations indicate that the non- 
circular motions produced by this instability provide the 
dominant mechanism for mass and angular momentum 
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transport in the disk (e.g. lKrumholz et alj|2007bll2009al) . 

Dozens of simulations of gravitational instability in 
disks have been publis hed, both for disks undergo- 
ing external accretion (fyorobvov fc Basul 120071 120081 
l2009t iKratter et a Hl201tt iMachida et al.ll2010l e.g .) and 

Ricd 12004. 2005 



for those in iso lation 
Kim &: Ostriker 



2007; ICai et al 



Lodato 

— 1 — 



20081: iCossins et 
20091; lAgertz et alJ 120091: see iDurisen et al.l 120071 for a 
review of earlier work). Based on these, several au- 
thors have presented one-dimensional time-dependent 
disk evolution models in which the effects of gravitational 
instability are approximated by an a prescription, with 
a obtained by fit s to simulation results or by general en- 
ergy arguments (IGoodmanl 120031; iHueso fc Guillo tll200 
Kratter et all 120081 ; IRice fc Armitagd 120091: IRice et al 
order-of-magnitude energy arguments 
to the case of galac t ic dis ks by 
Dekel et all (l2009bl). iKlessen fc Hennebelld (f2009h . and 
Elmegreen fc Burkertl (12010ft. In th e realrn of purely an- 



2010). Similar 
have been e xtende d 
(I2009bft. 



Li- 



alvtic work. lBertin fc Lodato! (11999ft present steady-state 
solutions for self-gravitating disk s with decayin g turbu- 
lence, while IRafikovl (|2009f ) and IClarkd pOOl derived 
steady-state accretion rates for disks in balance between 
radiative cooling and accretion-driven heating in proto- 
stcllar disks. 

While the analytic and one-dimensional models have 
provided a good understanding of the basic mechanism of 
gravitationally-driven turbulence and transport in disks, 
they also suffer from significant weaknesses. No an- 
alytic models published to date consider the case of 
gravitational instability-dominated disks that are time- 
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dependent rather than in steady state. Most previous 
work has been limited to a particular rotation curve (e.g. 
Keplerian or flat), to a disk of pure gas without stars, and 
to a disk that is vertically supported by thermal pressure 
rather than supersonic turbulence. As a result of these 
limitations it is not clear under what circumstances disks 
that are not in equilibrium can be expected to evolve into 
it, and it is not even clear what the equilibrium state is 
for a star-forming, supersonically turbulent disk such as 
a galactic disk. 

Our goal is to improve this situation by developing a 
first-principles theory for the evolution of a thin, super- 
sonically turbulent disk of star-forming gas in a state of 
marginal gravitational stability, driven by a specified rate 
of external accretion. The only assumptions we make are 
(1) that the disk maintains a state of marginal gravita- 
tional instability (Q — 1) at all times, and (2) the rate 
of energy loss due to radiative cooling can be parameter- 
ized as a certain fraction of the energy per crossing time 
of a disk scale height. We do not assume that the disk 
is in steady state or that it is characterized by any par- 
ticular rotation curve. From these assumptions, in Sec- 
tion [2] we derive equations describing the instantaneous, 
position-dependent rates of mass, energy, and angular 
momentum transport, and the time evolution of the disk 
surface density and velocity dispersion. In Section [3] we 
show that these equations admit an exact steady-state 
solution, and we derive the steady-state profiles of sur- 
face density, velocity dispersion, and transport of mass, 
energy, and angular momentum. In Section |4] we show 
that disks that are not in the steady state will evolve to- 
ward it, and that for high accretion rates this evolution 
occurs on an orbital timescale. Finally, in Section [5] we 
discuss the implications of our findings, and we summa- 
rize in Section [6] 

2. EVOLUTION EQUATIONS 

2.1. Mass, Angular Momentum, and Energy Transport 

We begin from the equations describing the evolution 
of a viscous fluid in a gravitational field that is in the pro- 
cess of turning its mass into collisionless stars. These are 
the equation of continuity, the Navier-Stokes equation, 
and the first law of thermodynamics: 
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Here p, v, e, and p, are the gas density, velocity, spe- 
cific internal energy, and pressure, p* is the rate per unit 
volume at which gas mass turns into stellar mass, ip is 
the gravitational potential, T is the viscous stress ten- 
sor, $ = T^(dvi/dxj) is the dissipation function, and T 
and A are the volumetric rates of energy gain and loss 
due to non-fluid (e.g. radiative or chemical) processes. 
Note that no terms associated with star formation ap- 
pear in the first law of thermodynamics or the Navier- 
Stokes equations because star formation does not alter 
the bulk velocity or specific internal energy of the gas. 

We consider a thin, axisymmetric disk centered on the 
origin lying in the plane z = 0. At every radius r the 



disk is characterized by a surface density £ and a total 
thermal plus non-thermal velocity dispersion a. The ma- 
terial orbits the origin with angular velocity and has 
a radial velocity ti r < Vf In Appendix lAl we show that 
for such a star-forming disk equations (JTJ, (jSJ, and ([3]) 
imply 
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is the viscous torque, j — rv^ is the specific angu- 
lar momentum, Q = v$/r is the angular velocity, ft = 
dlnvi/y/dhir, and Q = J T dz and C — J Adz are the 
vertically-integrated rates of non-fluid energy gain and 
loss. Equations Q, ([5]), and ([6]) are the standard equa- 
tions of mass, angular m omentum, and energy co nserva- 
tion for a thin disk (e.g. IBalbus fc Hawlevlll998D gener- 
alized to the case of a supersonically turbulent gas that 
is forming stars. 

If we assume that the disk is always close to radial 
force balance and that the potential varies slowly in time 
then we have dip/dr ss v\jr and dj/dt « 0. Using these 

conditions in equations (@|, ©, and ©, we arrive at the 
evolution equations for £ and a: 
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Finally, we note that the underlying assumption of our 
model is that we can represent the transport processes 
in a disk dominated by gravitational instability using a 
local viscosity. This assumption is only valid in certain 
circumstances, and this sets limits on the applicability of 
our model that we discuss in § 15.51 

2.2. Star Formation and Radiative Gain and Loss 

Equations ([5]) and © fully specify the time evolution 
of the system. The only physical approximations we have 
made thus far are that the disk is thin and axisymmet- 
ric, and that turbulent eddies on size scales of the disk 
scale height provide an effective pressure proportional to 
the square of the turbulent velocity dispersion. However, 
we have not yet determined the functions describing the 
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star formation rate £», the rates of radiative gain and 
loss Q and £, and the torque T. Since the physics in- 
volved in these terms is complex, we proceed in a simple 
parameterized manner. 

For the sta r formation rate, we note that both ob- 
servation (e g. IKrumholz fc Tanl l2 007l: IBigiel et aLllifflOl 
Evans et a.1.1 12009^) and theory (e.g. iKrumholz fc McKed 
2005t IKrumholz et al.l [20091 iPadoa n fc Nordlundl 120091: 
Hennebelle fc Chabrierll2'00l indicate that molecular gas 
forms stars at a rate of ~ 1% of its mass per free-fall time. 
We compute the free-fall time using the mid-plane den- 
sity; since the scale height is H = cr/fi, this is p — Efi/cr. 
Thus we adopt a star formation rate 



e ff SVG/O = £ff 



^G£ 3 ft 



1/2 



(10) 



where eg « 0.01. The true value of eg is likely to be 
slightly higher than this, because the clouds where stars 
form are generally denser than the mean midplane den- 
sity in the disk. Even with this correction, however, we 
can be confident that es < 0.1. Moreover, if a significant 
fraction of the ISM is atomic rather than molecular, the 
star formation rate is g reatly reduced (e.g. iWvder et al.l 
I2009t iBlanc et al.l 12009). in which case eg can be much 
smaller. More co mplex and accurate sta r formation laws 
are possible (e.g. IKrumholz et"aT1l2009bl ). but we will see 
below that such increased accuracy is not necessary at 
this stage. However, we do note that for a gaseous disk 
with Q — 1 (see Section HOI equation IT3"|) and a flat rota- 
tion curve (,5 = 0), equation (flU)) gives a star formation 
rate £* = 0.0067EJ7; in comparison, iKennicuttl ([1998) 
reports an observed star formation rate = 0.011SO, 
identical to within the errors Q 

For the energy loss rate, we note that numerous sim- 
ulations of turbulence show that it decays via radiative 
shocks on roughly a crossin g timescale (e.g. iStone et"al"1 
1998i lMac Low et al.l [l998l For the purpose of estimat- 
ing the loss rate we take the characteristic length scale of 
the turbulence to be comparable to the gas scale height 
H, so the crossing time is This is consistent with 

observations that show that turbulent power is domi- 
nated by the largest scales. We also limit our attention 
to galaxies where the total velocity dispersion is signifi- 
cantly in excess of the thermal value, since these are the 
only galaxies where one needs to explain the observed 
velocity dispersion by appealing to physics other than 
radiative balance. Thus we take the kinetic energy per 
unit area to be (3/2) Ecr 2 . The condition that the disk 
lose this amount of energy per crossing time of the disk 
scale height then reduces to 



(11) 



where rj is a dimensionless number of order unity. As 
a fiducial parameter we adopt r\ = 3/2, which corre- 
sponds to radiating away the full kinetic energy every 
scale height-crossing time. In disks where the velocity 
dispersion is primarily thermal rather than non-thermal, 
the loss rate will assume a different functional form (e.g. 



1 The coefficient reported in Kcnnicutt (1998) is 0.017 rather 
than 0.011. T he reduction to 0.011 comes fro m replacing the 
ISalpeteii {1953) IMF used in Kennicutt's work to a lChabriej l|2003 1 
IMF. 



Rafikov 2009) , but the remainder of our analysis will be 
unchanged. 

The gain rate is much more complex, since it involves 
turbulent motions generated by star formation. Since 
we are interested in systems where gravitationally-driven 
turbulence dominates, however, we make the extreme as- 
sumption that Q = 0. We return to the question of the 
real value of Q in Section 15.31 



2.3. Gravitational Stability and the Torque Equation 

We now turn to the central hypothesis of our model, 
which is that the a self-gravitating disk will adjust its 
torque, and therefore its radial mass flux, so as to remain 
in a state of marginal stability. T his hypothesis has als o 
been investigated in the models of iBurkert et al.l (|2009j) . 
For a disk of g as plus stars, t he parameter that describes 
its stability is (|Rafikovll200ll ) 



Q(q)- 1 =2Q 



-il 



where 



'Ml 2 ) 



q 2 R 2 ' 
(12) 



R=— , (13) 



and Iq is the Bessel function of order zero. Here £* 
and <7* are the surface density and velocity dispersion 
of stars, k = [2(/3 + l)] 1 / 2 ^ is the epicyclic frequency, 
and q = fcer*/ 'k is the dimensionless wavenumber of the 
mode in question. Modes for which Q(q) < 1 are un- 
stable. Note that iRomeo et al.l (|2010l ) have proposed a 
generalization of this condition for the case of gas with a 
scale-dependent velocity dispersion, as expected for tur- 
bulence, but for simplicity we use the iRafikovl (|2001h cri- 
terion instead. 

It is not generally possible to find the minimum value of 
Q(q) analytically. However, we can obtain a significant 
simplification if we focus on the most interesting cases 
for gravitationally-driven turbulence. In disks at high 
rcdshift there has not been time for the stars and gas to 
evolve so that their velocity dispersions are very different 
- see § 15.21 for a further discussion of this point. For 
these disks we therefore adopt er* = a. In this case the 
expression 



7rG(£ + £*) 



(14) 



is accurate to better than 7%. With this approximation, 
the condition for a disk to remain marginally stable at 
Q = 1 becomes 



~~fo~dt + dZ~dt 

_lda_ 1_ 

~~u ~0t ~ £ + 



dQ dT, t 
f <9£* dt 
<9£ <9£* 
~dt + ~dt 



(15) 
(16) 



Plugging in our expressions for the temporal derivatives 
of (j, E, and £*, we obtain an equation that describes 
the torque required to maintain Q = 1: 



dr 2 dr 
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with 
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where / 9 = E/(E + E„) is the gas fraction in the disk, 
and we have used Q = 1 to replace any dependence on 
E with a dependence on a and / 9 . We refer to this as 
the torque equation. Note that in deriving equation (JTTJ) 
we have implicitly assumed that stars do not migrate 
radially from their formation locations. If we were to 
make the opposite assumption, that stars and gas move 
together, then the combined gas plus star disk would 
act essentially like a purely gaseous one, except that 
the stars would be dissipation free. The corresponding 
torque equation is simply equation (|17j) with r\ replaced 
by ijfg, and f g — > 1 in all other terms. 

The other interesting location to consider for 
gravitationally-driven turbulence is in outer parts of 
present-day disks. In these regions star formation occurs 
at a negligible rate, and a* S> o\ In this case Q ss Q g , 
and a little calculation shows that the torque equation 
is simply (JTTJ) with f g — 1 everywhere. The stars sim- 
ply become irrelevant for gravitational stability. Thus 
equation (1171) . with the appropriate choice of f g and 77, is 
capable of representing both a present-day outer galaxy 
disk and a high redshift disk. 

To finish specifying the torque, we must provide 
boundary conditions for equation (fTTj) . One boundary 
condition comes from fixing the external accretion rate 
onto the galaxy to a value M = M cxt (so dT/dr = 
—v^{8 + l)M ex _t) at the disk outer edge, r = R. We imag- 
ine the inflow rate at this radius to be set by cosmological 
infall, which occurs at a rate unaffected by what happens 
in the disk. The other boundary condition depends on 
how we handle the inner boundary. If we imagine that 
our model applies all the way to r = 0, we must find a 
solution that remains regular in the vicinity of the sin- 
gular point there. We show below that there does exist 
a steady-state solution that is regular at r = and has 
the specified accretion rate M = M cxt at r = R. Outside 
of that steady state it is not possible to simultaneously 
have regularity at r = and an externally-imposed ac- 
cretion rate at r = R, so we must instead truncate the 
model at some radius tq > 0, where we imagine that a 
stellar bulge or some other non-disky structure forms. In 
that case we require that the torque have a small value 
at r — ro, so that the inner bulge region does not do 
work on the disk. 

The evolution of the system is now fully specified. At 
any given time, equation (1171) specifies the viscous torque. 
That torque in turn sets the time evolution of E, cr, and 



E* via equations ©, and (jTUJ) . Before moving on, 
however, we pause to point out some of the important 
physical properties embodied in equation (|17[) . First, 
the star formation rate does not appear explicitly in the 
torque equation. This is because for cr* ss a changing 
gas into stars does not significantly affect the stability of 
the disk, and with our approximate form for Q it does 
not affect the stability at all. Star formation enters the 
problem solely through its effects on f g . Second, if 77 = 
then clearly T — is a solution for Af oxt = 0. Physically, 
this represents the fact that, if there is neither dissipation 
of turbulence nor accretion, then the disk can remain 
marginally stable without any mass transport. 

2.4. N on- Dimensional Equations 

It is helpful at this point to non-dimensionalize our 
equations and derive some characteristic numbers. If we 
make a change of variables r = xR, a = sv^(R), = 
uvj,(R), and T = rM oyi tV ( f > (R)R, then we can rewrite 
equation (JTTJ) as 



t" + hxr' + h a r = H 
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h ■ 



hi 



H 
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(22) 

(23) 
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subject to the boundary condition t' = —B — 1 at x = 1. 
Here primes indicate differentiation with respect to x, 
and we have defined 



X 



GMoxt 

MrJ 3 



(26) 



The form of equation (|22|) is instructive. The coefficients 
ho and hi appearing on the left-hand side depend on the 
current state of the disk without reference to external 
accretion or turbulent dissipation. Those affect only the 
inhomogeneous term H on the right-hand side, which is 
proportional to r\jx- We may view H as the driving term 
for the system, with more rapid dissipation of turbulence 
(i.e. larger 77) and lower accretion rates (lower x) both 
tending to produce larger torques. 

We can similarly non-dimensionalize the evolution 
equation for the gas surface density and velocity dis- 
persion by defining E = SM C xt/(v<p(R)R) and t — 
T[2ttR/v c ( > (R)], so that the evolution equations become 
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we have defined £* = S^M^t/ {v^RjR) in analogy with 
S, and the Q — 1 condition in dimensionless form is 



Clearly 



y/2(P + l)us 
nXx(S + S*) 



1. 



(30) 



In these units the orbital period is 1, and the fraction of 
the disk mass that accretion adds per orbit is (S + S*)^ 1 , 
where the angle brackets indicate an average over the disk 
area. 

3. STEADY-STATE DISKS 
3.1. The Steady-State Solution 

Having derived the basic equations that govern the sys- 
tem, we now search for steady-state solutions, which we 
can obtain analytically. A true steady state is obviously 
not possible in a real disk that undergoes mass accretion 
and star formation, but we can find solutions which are 
steady for time periods that are short compared to the 
star formation and accretion timescales. We will there- 
fore set eff = 0, and in Section |3~21 we will check that this 
is a reasonable approximation. In this case a steady state 
solution is one for which dM/dr — 0, or in dimensionless 
form 

d 



dx 



{0 + l)u 



= 0. 



(31) 



Combined with the boundary condition that t' = — (0 + 
1) at x = 1, this implies that the steady state solution 
is t' — —{fi + l)u. One can immediately verify using 
equation ([27)1 that such a torque gives dS/dT = 0. 

To make further progress toward an analytic solution, 
we concentrate our attention on cases where has a con- 
stant value. This is a reasonable limitation, since most 
galaxies have flat rotation curves (0 = 0), while Kep- 
lerian disks have = —1/2. For constant 0, we have 
u = x@ , and we can analytically integrate the steady 
solution for t' to obtain r = — x l3+1 + c, where c is a 
constant to be determined by the regularity condition 
at the inner boundary. In Appendix [B] we provide this 
analysis for two of the most physically important cases, 
= (flat rotation) and = —1/2 (Keplerian rotation), 
which shows that c = in both those cases. Of course 
cannot have a constant value of or —1/2 all the way 
to r = 0, since the rotation velocity would diverge, but 
our approximation is appropriate in cases where has a 
constant value over a large dynamic range in radius, and 
changes only at small radii where a bulge forms (in the 
case of a flat rotation-curve galaxy) or a boundary layer 
joins a disk to a star (for a Keplerian disk). 

Given the value that r must have in order to produce 
a steady state, we can plug it into the torque equation 
to determine the corresponding disk properties that are 
required. For a flat rotation curve, = 0, we have 

1 



(3/ s - I)* 2 * 2 
( 2 3/2 / 9 



I ?• (32) 



\3/g - 1/ \XJ X 

and the steady-state condition r = — x then implies 

j = 2y/2f g ris 3 - x 
bsxx 



(33) 
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is an exact solution, and numerical integration shows 
that all solutions converge to this value very quickly at 
x < 1 regardless of the value of s at x = 1. The corre- 
sponding surface density and inward velocity of the ma- 
terial are, from equations (|A2I) and (fT4"|) . 
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and the corresponding d i mensi onless viscosity param- 
eter iShakura fc Sunvaevl ([1973) and viscous evolution 
timescales are 



GM 2V2 
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v r(R) \VX 



Mj\ — , (38) 



where t OI \, = 2-kR/v^R) is the outer disk orbital period. 
(We omit a factor of Q in the equation for a because 
it is set to unity in our model.) For the corresponding 
Keplerian case (/? = — 1/2) it is easy to verify by a similar 
procedure, and using the analysis of the singular point 
provided in Appendix [Bl that the steady solution is r = 

-VS. s = [3x/(4t7/ 3 )] 1/3 . 

Our result is easy to understand intuitively. For a 3> 
c s , the rate per unit mass at which the turbulence decays 
is 77<7 2 f2, so the turbulent decay time is of order the or- 
bital timescale times r\. To keep the velocity dispersion 
constant, mass must move inward at a rate such that the 
decrease in gravitational potential energy balances bal- 
ances this radiative loss. For a flat rotation curve, the 
decrease in potential involved in moving from radius ro 
to radius r is ln(r/ro), so inward drift at a velocity v r 
causes the potential energy per unit mass to decrease at 
a rate v^(v r /r). Equating the rates of dissipation and 

energy increase gives i]a 2 fl = vi(v r /r), and it follows 

immediately that v r = rja 2 jv$. 

It is also worth noting that this result is very similar to 
that of iGamml e (2001), who finds that, in steady state 
in a disk that cools on a timescale r c , gravitationally- 
induced turbulence produces an effective viscosity a = 
[7(7 — l)(9/4)f2r c ] -1 . Our effective "cooling time" for 
the supersonic turbulence is t c — 1/(77^), and the re- 
mainder of our result differs from his only in that we 
have modeled disks with a stellar component, and that 
our pressure and internal energy are appropriate for su- 
personically turbulent motion on scales comparable to 
the disk scale height, rather than for an adiabatic gas 
described by a polytropic equation of state. The former 
introduces a dependence on f g , and the latte r produce s 
the slight change in leading coefficient. That Ga mmiel 's 
result can be obtained on energetic grounds, rather than 
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by computing stresses as i n his d eriva tion, has also b een 
pointed out bv lRice et all (120051 ) and lRafikovl (|2009T >. 

3.2. Conditions for Steady-State 

In deriving the steady solution we have ignored the 
change in total disk mass due to accretion. More subtly, 
our steady solution has constant M all the way in to 
r = 0, so mass effectively vanishes through the origin. 
These approximations are only reasonable for time scales 
over which the total disk mass changes little. We can 
define the accretion timescale for our steady solution as 

r R 27rr(E + E„)rfr _ t olb _ t olh 

Mext " *(vf g X 2 ) 1/3 ~ 2n V f g sf 

(39) 

Note that t acc is the time required for external accretion 
to double the disk mass, and is distinct from the viscous 
accretion time defined by equation (|38|l . To avoid confu- 
sion we will always refer to that quantity as the viscous 
timescale, and use the accretion time solely to refer to 
the mass-doubling time produced by external infall. For 
our fiducial rj = 3/2, we have 

t acc ~ 1.3/f 1 /3 x -a/3 torbj (40) 

where xo.i = x/0.1. Similarly, we have neglected star 
formation, which is only reasonable for time short com- 
pared to the star formation timescale. We can define 
the star formation timescale as the mean ratio of star 
formation rate to gas surface density in the disk: 

tsF = — / (27rx) S dx = tmh . — . 

7T J J dS,/dT (1627r 2 ) 1 /4/i/2 eff 

(41) 

For our fiducial eg — 0.01, this gives 

i SF ~ 16f- 1/2 t OTb . (42) 

In comparison, the characteristic evolution timescale 
for the disk should be the viscous time defined by equa- 
tion ([38]). since this is the time required to drain the 
material in the disk and replace it with newly- accreted 
material. For our fiducial 77 = 3/2 and a flat rotation 
curve, this is 

i visc = 1.3/»/ 3 Xo.? /8 *orb. (43) 

Comparing the viscous time to the timescales relevant 
for accretion and star formation, we have 

j^=f 3 (44) 

«ac.c 

^=(5) ,,4 |^=»-0<V.r, (45) 

for our fiducial values of 77 and eg. Our neglect of 
accretion-induced changes in the disk mass and star 
formation-induced changes in the gas fraction is reason- 
able when these two ratios are < 1. Clearly the require- 
ment that iviscAacc J$ 1 is always satisfied, although per- 
haps only marginally if f g is large. The requirement that 
iviscAsF < 1 is satisfied as long as x > 10 -3 , or longer 
if the disk is non-star- forming (eg = 0). Thus we expect 
that our steady state model is reasonable under these 
conditions. 



4. TIME-DEPENDENT SOLUTIONS 

We now consider a disk that is initially out of steady 
state, with a specified initial value of \ an d fg- For 
the reasons discussed in the previous section, we take 
eg = and x an d f g as constant. We start each calcula- 
tion from an initial velocity dispersion specified on a grid 
of N x cells, logarithmically spaced. The grid runs from 
x = xq to 1, where xq is close to zero. For the bound- 
ary conditions, we cannot simultaneously require that r 
obey the regularity condition derived in Appendix [B] at 
the inner boundary and that t' = — (/9+ l)u at the outer 
boundary; this would amount to applying three bound- 
ary conditions to a second-order ODE, and for general 
choices of s no solution for r exists that satisfies all three 
conditions. Instead, we continue to fix r' (and thus the 
accretion rate) at the outer boundary, and at the inner 
boundary we require that the torque be r = — xq. This 
choice amounts to requiring that the work done by gas 
in the region x < xq on the computational domain at 
x > xq approaches as io — » 0, while the mass flux 
through the inner boundary at x — xq is allowed to vary 
freely. This is a good approximation to a disk being fed 
from the outside a fixed rate and that has an inner bulge 
region or an inner boundary layer that is stress free, but 
which can accept mass at varying rates. Note that we do 
not set r = at x — xq exactly, because this is inconsis- 
tent with the steady state solution. 

With this setup, we evolve the system accordin g to 
equation (j2"8l using the algorithm given in Appendix QB- 
Note that we find that the unmodified evolution equa- 
tion (|28l) for s is numerically unstable to the growth of 
small oscillations on the grid scale. We damp these by 
adding a small amount of viscosity to the disk evolution, 
implemented in a manner that maintains exact energy 
conservation. 

In Figure [1] we show the evolution of disks with (3 = 
(flat rotation curve), f g = 0.5, r\ = 3/2, xo = 0.1 and 
N x = 500. The left panels show runs with \ — 0.1, and 
the right panels show runs with \ = 0.01. The top row 
shows disks with an initial velocity dispersion sq equal 
to either twice the equilibrium value s oq given by equa- 
tion (|34l) , and the bottom row shows disks with an initial 
velocity dispersion that is half this value. As the figure 
shows, all disks evolve toward the equilibrium solution 
found in Section [3] very rapidly, regardless of whether 
they start with initial velocity dispersions smaller or 
larger than s oq . The runs with % = 0.1 halve their dis- 
tance from the equilibrium solution in less than a quarter 
of an outer orbital period, and converge to within 10% 
of the equilibrium solution within one full outer orbit, 
after which point they are essentially static0 in fact, the 
time required to reach equilibrium seems to be a factor 
of ~ 2 — 3 less than our naive estimate that it should 
the viscous time t v i sc . This may be because 7j v isc is an 
estimate of the time for the material to reach zero ra- 
dius, while in our case the material need not travel as far 

2 A Mathcmatica program to implement this algorithm is avail- 
able at http://www.ucolick.org/~krumholz/downloads.html 

3 The x = 0-1 cases miss the equilibrium value very slightly 
and converge to a velocity dispersion that is a few percent above 
it. This is an artifact of the small viscosity we require in order to 
maintain numerical stability in this run. The \ = 0.01 cases are 
stable with a somewhat smaller viscosity, so the deviation from the 
exact equilibrium is unnoticeable for them. 
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Fig. 1. — Time evolution of the velocity dispersion s as a function of radius x for runs with \ = 0.1 (left column) and \ = 0.01 (right 
column), and with initial velocity dispersions sq = 2s eq (top row) and so = ^eq/2 (bottom row) at all radii. Here s oq = Ix/ivfg)] 1 / 3 /v / 2 
(equation 134 \ is the analytically-computed equilibrium velocity dispersion, and these runs use rj = 3/2, /„ = 1/2. In each plot, the black 
curve labeled r = is the initial velocity dispersion, the red curve is the velocity dispersion at r = 0.25 (i.e. after 0.25 outer orbits), and 
each subsequent curve after that represents a factor of 2 increase in time: r = 0.5 (green), r = 1 (blue), r = 2 (purple), and, for the 
X = 0.01 runs, r = 4 (aqua). The thick dashed black line is s = s aq . 



to set up an equilibrium velocity dispersion and surface 
density profile. 

Thus the time required to reach equilibrium is well be- 
low the accretion time i a cc = 1.6t or b (for f g — 1/2) we 
computed in Equation (j40|) . and is far less than the star 
formation time isF = 20i or b given by Equation (142 p . in 
accord with our analytic expectations. We therefore con- 
clude that disks with % = 0.1 converge to the equilibrium 
configuration on a timescale short compared to both the 
accretion and star formation timescales. For \ = 0-01 
the convergence to equilibrium is slightly slower, but the 
runs are within 10% of equilibrium by 2 outer orbits, and 
are within ~ 1% of equilibrium by 4 outer orbits. Since 
the accretion timescale is 7.6 orbital times for x = 0.01, 
and the star formation time is 20i or b, these runs too reach 
equilibrium fast compared to i acc or tsv- 

Thus we have demonstrated that the time-independent 
solution we obtained in Section [3] is not only an exact 
equilibrium, it is an attractor toward which initially out- 
of-equilbrium disks will converge. As long as \ > 10 -3 
(or for arbitrarily small x m non-star- forming disks) , this 
convergence occurs on a timescale short compared to ei- 
ther the either the star formation timescale or the accre- 
tion timescale over which the disk mass changes appre- 
ciably. This means that our earlier decision to neglect 
both star formation and changes in the rotation curve 
due to accretion is reasonable, and that a disk that is 
out of steady state will converge to its time-independent 
equilibrium state faster than either star formation or ac- 
cretion can alter that equilibrium. It also implies a vast 
simplification in comparing to observations: since con- 
vergence to equilibrium is fast, we can generally assume 
that the velocity dispersions and surface densities of ob- 




FlG. 2. — Disk velocity dispersion a versus redshift z for halos of 
mass Mh 12 = 0.1, 0.3, and 1.0 (blue, red, and black lines) and gas 
fraction f g = 1/2 or 1 (dashed and solid lines). The plot uses our 
fiducial rj = 3/2 and assumes /), q.is = 1/2, i.e. half the infalling 
baryons arc gas and half are stars. The f g = 1 case is appropriate 
for systems where there no stars or where cr* S> o", such as present- 
day galactic disks, while the case f g = 1/2 is appropriate for high- 
redshift disks where tr* Rj a. 

served disks reflect the instantaneous equilibrium state 
dictated by their current external accretion rates and gas 
fractions. Of course we still have not included star for- 
mation feedback, a topic we approach in Section 15.31 

5. DISCUSSION 

5.1. Cosmological Evolution of the Velocity Dispersion 
of Galactic Disks 

We have now shown that disks dominated by 
gravitationally-driven turbulence rapidly converge to an 
equilibrium state in which their velocity dispersions are 
determined by their gas fractions and accretion rates. 



Since this convergence happens on an orbital timescale, 
most galactic disks should be found near their equilib- 
rium state. We can use this result, coupled to a simple 
model for how galaxy halos accrete mass, to study the 
evolution of disk velocity dispersions over cosmic time. 
Using Press-Sche c hter f its to dark matter simulations, 
iNeistein fc Dekel ([20081 ) estimate the mean dark mat- 
ter accre tion rate onto halos of a given mass at a given 
rcdshift. iBouche et al.l ([2009D extend this to give an esti- 
mate of the gas accretion rate onto the disk at the center 
of the halo, which we adopt: 

M g = 7.0^/6,0.18^2 (1 + zf- 2 M yr-\ (46) 

where z is the redshift, /^o.is is the gas mass fraction of 
the infall divided by 0.18, the universal baryon fraction, 
Mh,i2 is the halo mass in units of 10 12 M©, and ei n is the 
fraction of gas entering the halo that reaches the galactic 
disk rather than being shock-heated and joining the halo. 
This is approximately given by 

_/0.7/(z), M M2 <1.5 U7] 
in ~ \ 0, M hA2 > 1.5 ' 

where f(z) is a function that is linear in time and varies 
from unity at z — 2.2 to 0.5 at z — 10 Inserting M g from 

equation (j46]l for M ext in equation (j34|). we are able to 
evaluate the expected velocity dispersion of gravitational 
instability-dominated galactic disks as a function of halo 
mass and redshift. We do so in Figure [2j 

Examining the plot, we see that for a Milky Way- like 
halo (M h ,i2 = 1, Kue et al.1 120081 ) where ct* > a (so 
that the f g = l case applies), we predict a typical veloc- 
ity dispersion of 10.7 km s _1 . While this is very slightly 
higher than the value of a ~ 8 km s -1 observed in typical 
Milky Way-like disks today (jBlitz fc Rosolowskvl 12004" 
iDib et all 120061 and references therein), the agreement 
is quite good given our purely analytic model. Our re- 
sults are quite similar t o the numerical ones obt ained by 
iKim fc Ostrikerl (12007ft and lAgertz et al.l (|2009| ). More- 
over. as lAgertz et al.l point out, gravitationally-driven 
turbulence has the advantage that it can operate even 
in the outer H 1 disk where there is very little star 
formation, so mechanisms such as supernovae that are 
invoked to explain turbulence in the inner disk (e.g. 
Ide Avillez fc Breitschwe"rdtll2007t Uoung et al.l I2009Q are 
unavailable. We also predict lower velocity dispersion 
in smaller halos, and this appears to be consistent with 
the somewhat lower H 1 velocity dispersions seen i n 
dwarf galaxies ([Walter et al.l 120081: iChung et all 120091 ) . 
Finally, however, we do note that there are alterna- 
tive models to explain outer disk turbulence, includin g 
magnetorotational instabi lity lSellwood fc Balbusl (|1999f) : 
Piontek fc Ostrikerl (12007ft and accretion of clumpy gas 
Sant illan et al l ([20071 ). 

Within the same framework we are able to explain the 
large velocity dispersions of 20 — 80 km s " 1 found in 
galac tic disks found at redshifts ~ 1.5 — 3 ([Cresci et al.l 
2009). The observed galaxies likely correspond to ~ 10 12 
Mq halos. For redshifts in this ra nge and f„ ~ 1/2 , 
typical of galaxies at that redshift (Daddi e t al.l 120091 : 

4 We compute the time as a function of redshift, and all other 
cosmology-dependent quantities, using Q m = 0.28, Qa = 0.72, and 
h = 0.70. 
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Fig. 3. — Ratio of disk maximum circular velocity V max to ve- 
locity dispersion a as a function of redshift z for halos of mass 
Mh 12 = 0.1, 0.3, and 1.0 (blue, red, and black lines, top to bot- 
tom) and gas fraction f g = 0.1, 0.5, or 1 (dot-dashed, dashed, and 
solid lines). All parameters are the same as for Figure [2] 

iTacconi et al.l I2010D we predict typical velocity disper- 
sions of 30 — 50 km s _1 , with fluctuations at the factor 
of ~ 1.5 level, corresponding to the expected factor of 
~ 3 level variations in the accretion rates of halos at the 
same mass and redshift. This is in good agreement with 
the observations. 

It is also instructive to compare the velocity dispersions 
we predict to the expected rotation velocities of galactic 
disks. We compute the approximate virial velocity of a 
halo as a function of mass and redshi ft following the ap- 
proxim ation given in Appendix A2 of iDekel fc Birnboiml 
(pPM ). and we take the maximum circular velocity to 
be 1.2 times this bas ed on fitting the ze ro point of the 
Tully-Fisher relation ([Dutton et al.ll2007ft . With this ap- 
proximation, we plot Vmax/cr in Figure [3l We see that 
accretion-driven turbulence naturally produces the tran- 
sition from disks with V max /cr ~ 5 found at redshifts > 2 
to disks with V maj! _/cr ^ 20 — 25 found today. 

Interestingly, we find that there is little dependence of 
Knax/c on halo mass. Instead, the primary dependence 
is in f g , the gas mass fraction; analytically, V max /cr oc 

— 1/3 

f g . Thus the most gas-dominated systems (or old 
galaxies that have 3> <r) have the largest V max /cr, 
while gas-poor systems have smaller V ma ^/a. This sug- 
gests that the range of Vm^/a seen for galaxies at z ~ 2 
by the SINS survey ([Cresci et al.ll2009T ) represents a se- 
quence in gas fraction. The dispersion-dominated galax- 
ies should on average be comparatively gas poor, while 
rotation-dominated ones should be gas rich. Of course 
fluctuations in accretion rate can also cause changes in a, 
so detecting this effect will require samples large enough 
for this noise source to be averaged out. Nonetheless, it 
seems likely that data to test this prediction will become 
available in the next few years. 

5.2. High- Redshift Galaxies 

It is particularly interesting to apply our models to the 
z ~ 2 — 3 galaxies ob served by lElmegreen et al.l (12004 
2005), [Genzel et all (12008ft . iForster Schreiber et all 
(|2009t ). ICresci et al.l (|2009t ). and others, since these 
are thought to be examples of strongly gravitational 
instability-dominated disks. We first note that, in the 
rcdshift range z — 2 — 3 for halos of mass Mh = 10 12 Mq, 
thought to be typical of the observed systems, the models 
shown in Figures [2] and [3] give x = 8x 10~ 3 — 1.1 x 10~ 2 . 
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For these values of x an d g as fractions f g = 1/2, using 
equations (|40|) and (1421 the ratio of star formation time 
to accretion time ts^/t acc = 1.8— 2.6, so the star forma- 
tion rate is roughly 1/2 — 1/3 of the total accretion rate. 
Given the uncertainties in this model and the dispersion 
in expected accretion rates, for simplicity we can simply 
adopt M» « M cxt - Since the star formation rates are 
observed (and have typical values ~ 100 M Q yr _1 ), we 
can plug them into our model in place of M ext in order 
to predict disk properties. 

Doing so, we find that the redshift 2 — 3 disks should 
have velocity dispersions (from equation [ 



47 km s- 1 f- 1/3 Ml{ 3 



100' 



(48) 



where M* t too = M*/100 Mq yr _1 , independent of their 
maximum rotation velocities V max . Thus galaxies of sim- 
ilar star formation rate and gas fraction should have 
the same a independent of V max . The viscous accretion 
timescale required for the gas at the edge of one of these 
disks to reach the center is (from equation |3"5|) 

t visc « 600 Myr fi /3 Ri Vj£M;$, (49) 

where Rio = R/10 kpc and V200 = V max /200 km s _1 , 
and the gas mass is (from equation [33]) 



M g = Rv 4 



1/3 



= 3 x 10 10 M Q f g '*Rx Q V 2aQ Ml%. (50) 

The ratio of baryonic to dynamical mass within the disk 
region R is 



M dv 



3 f -1 / 3 !/ -1 M 1 / 2 

u -°Jg K 2Q0 iW *,100- 



(51) 



To the extent that these quantities have been observed, 
they are in good agreement with the results of our model. 

It is also useful to verify that our approximation that 
<t» sa a is valid for these galaxies. Once stars form, tran- 
sient spiral structures will dynamically heat them until 
the ste llar disk becomes stable aga i nst further spiral pat- 
terns (ISellwood fc Carlberd [1981 ICarlberg fc Sellwoodl 
1985). The characteristic timescale for this heating is 
[(Qiim — Q*)/T~)t or h, where Qu m « 2 is the limiting value 
at which the disk becomes stable against spiral pertur- 
bations, is the current Toomre Q parameter for the 
stars, and numerical experiments show that r ~ 4 — 5 
for i or b evaluated at one disk scale length. If we assume 
that the stars are born at Q* = 1, then the characteristic 
timescale over which the heat is 



750 Myr R 10 V^, 



(52) 



where we have taken the radial scale length to be half of 
R. In contrast, the time required to double the stellar 



1 - fg ( M < 



fg 



) = 300 Myr 



,1/3 
J 9 



RiqVwqM, 



-2/3 



(53) 



Thus we see that the time required for stars to increase 
their velocity dispersion via spiral structure is generally 



comparable to or longer than the time required for a 
new generation of stars to form with the same velocity 
dispersion as the gas. Our approximation that the stellar 
population has the same velocity dispersion as the gas in 
these galaxies is therefore reasonable. 

5.3. Effects of Stellar Feedback 

In our idealized models, we have neglected the influ- 
ence of stellar feedback by setting Q = 0. This is ob- 
viously reasonable if we are concerned with the outer 
parts of a galactic disk where there is no star forma- 
tion, or a protostellar disk where at most a few stars will 
form. It is not reasonable for the centers of present-day 
galactic disks, where supernovae are clearly important. 
It is questionable whether stellar feedback is important 
in ULIRGs or the high surface density galaxies found 
in the early univer se. Supernovae are no t effective in 
such environments ([Thompson et al.l l2005t Uoung et al.l 
2009), but stellar radiation pressure may be. Whether 
it can actually drive the observed velocity dispersions in 
high redshift galaxies is a matter o f debate ([Murray et al.l 
2009; IKrumholz fc Matznerl[2009h . 

In those situations where feedback is significant, we 
can qualitatively see how it would change our results by 
noting that adding a non-zero Q to our equations would 
have roughly the same effect as lowering r\. Physically, 
if star formation injects turbulence into the ISM at an 
appreciable rate, this is equivalent to reducing the rate 
at which turbulence decays - we effectively increase the 
"cooling time" of the disk. Consulting equations (|34[) 

([3"7| . we see that the effect of this is to increase the 
velocity dispersion and surface density in the equilibrium 
state, while reducing the radial velocity and the rate of 
angular momentum transport. 

We caution that this analysis is only valid as long as 
the feedback is not too strong. In particular, we require 
that C remains larger than Q for a Q = 1 disk, and that 
the turbulent stresses created by the feedback mecha- 
nism are significantly weaker than the stresses induced 
by gravitational instability-driven turbulence. If the first 
requirement is not met, then feedback will drive the ve- 
locity dispersion up to the point where Q > 1, and the 
gravitational instability will shut off. If the latter condi- 
tion fails, then gravitational instability will continue, but 
our calculation of the transport rate will not be correct 
because we have not included stresses induced by feed- 
back. Even if both requirements are met, our analysis of 
feedback effects should be regarded as qualitative rather 
than quantitative. Energy injection Q appears on the 
right hand side of the torque equation with the opposite 
sign as C, but their functional dependence on other disk 
parameters (surface density, velocity dispersion, etc.) are 
almost certainly different. The exact effects of feedback 
will depend on how energy injection varies with these 
quantities, which will in turn depend on the type of feed- 
back and the physics of the ISM. 

5.4. Protostellar Disks 

Although we have focused our discussion thus far on 
galactic disks, our model applies for arbitrary rotation 
curves, gas fractions, and infall rates, so we can apply 
it equally well to protostellar disks. To understand the 
expected levels of turbulence in protostellar disks, we use 
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the p arameterization of infall due to lKratter et al.l ((20081 
2010), who introduce the dimensionless numbers: 



£ 



GM V - 



r = 



(54) 



where c s .d is the sound speed in the disk, M*d is the 
total mass of the disk and star and f2& in is the Keple- 
rian orbital period of the infalling material. Physically, £ 
represents the ratio of the external accretion rate to the 
maximum rate (~ c 3 d /G) at which a stable disk can pro- 
cess material, while T represents (up to a factor of 2tt) 
the fraction by which the disk plus star mass changes 
per outer disk orbit. Indeed, since v^(R) = R£lk,in = 
y/GM+d/R, with a little algebra it is easy to show that, 
in the case of a Kepleria n disk consisti ng entirely of gas, 
our x simply reduces to lKratter et all s V parameter. 

With this understanding, we can explain the observa- 
tion by iKratter et al.l (|2010l ) that, in their simulations, 
the typical velocity dispersion of disks that do not frag- 
ment is comparable to the disk thermal sound speed (see 
their Figure 8). For a purely gaseous Keplerian disk, our 
model gives s = [3x/(4r?)] 1 / 3 , and IKratter et all show 
that the disk sound speed is related to the Keplerian ve- 
locity at the disk edge by c s ,d/v^{R) = (r/^) 1 / 3 (their 
Equation 18). Combining these two results, the expected 
Mach number of the accretion-driven turbulence is 



M = 



C s ,d 



Cs,d 



1/3 



(55) 



Since fragmentation is avoided only for disks with £ of no 
more than a few, we can take £ ~ 1, and it immediately 
follows that the expected Mach number M. ~ 1. 

We can apply a similar analysis to real protostellar 
disks: the Mach number of the turbulence in these disks 
should follow equation (f55|) . This means that disks ac- 
creting with £ ~ 1, corresponding to M cxt ~ 10~ 5 
Mq yr _1 for typical out disk temperatures T ~ 50 K, 
should have disks whose turbulent velocity dispersions 
are roughly transonic. This state should prevail dur- 
ing the majority of the main accretion phase. Once the 
main accretion phase ends and the accretion rate drops, 
the turbulent velocity dispersion should drop to subsonic 
values. This represents another prediction from our anal- 
ysis: class and class I protostars should have disks with 
transsonic turbulent velocity dispersions, while class II 
and class III sources should have subsonic turbulent ve- 
locity dispersions. As ALMA comes online in the next 
few years and provides resolved molecular line maps of 
prot ostellar disks at a variet y of stages in their evolution 
(e.g. iKrumholz et~a l. 2007a), we will be able to test this 
prediction. 

5.5. On the Validity of a Local Viscous Approximation 
for Gravitational Instability- Induced Transport 

The central approximation we make in our model is 
that transport of mass, angular momentum, and en- 
ergy produced by gravitationally-driven turbulence can 
be represented with a local viscous stress tensor. The 
validity of this approximation ha s been the subject of 
great debate in the past decade. iBalbus fc Papaloizoul 
(|1999ft show that self-gravitating disks cannot in general 



be modeled with a viscous formalism, but that such an 
approximation may be reasonable for disks near Q = 1, 
the condition that we adopt throughout this work, and 
that appears to apply to the galactic and protostellar 
disks we are interested in studying. Based on a com- 
bination of an alytic arguments and local simulations, 
IGammid (pOOlh argues that a local prescription is appli- 
cable to Q = 1 disks that are suffic iently thin, s < 0.12 , 
and more recent global simulations (lLodato fc Ricell2004 
[2001 IBolev et all IpOOl ICossins et al.l 120091 ) generally 
support this result. IGammier s condition for a local trans- 
port approximation to apply is well-satisfied for galactic 
disks at redshifts z < 2 (§ 15.11) and for non-fragmenting 
protostellar disks (§ I5.4p . It is marginally violated for 
the observed disks at z ~ 2 (§ 15.21) . suggesting that our 
model should be considered with some caution for them. 
At a minimum the thickness of these disks likely pro- 
duces different fragmentation be havior than a standard 
thin disk analysis would suggest (jBegelman fc Shlosmanl 
l2009h , 

6. SUMMARY 

In this paper we derive the basic evolution equations 
for a disk of gas and stars kept in a state of marginal 
gravitational instability by a combination of external ac- 
cretion, inward migration of gas, and decay of turbulent 
motions due to radiative shocks. In such a disk, we use 
the equations of conservation of mass, angular momen- 
tum, and energy to derive an equation (|22p that char- 
acterizes the instantaneous rates of mass and angular 
momentum transport required to maintain the state of 
marginal stability, and we show that this equation has an 
analytic steady-state solution in which the disk velocity 
dispersion (Equation IMl) . surface density (Equation I55"j) . 
and rates of transport (Equations l5B1andl57]) through the 
disk are determined by the rate of external infall onto the 
disk and the gas mass fraction within it. We show that 
disks converge to this steady state on timescales of order 
the orbital time, much less than the time over which ei- 
ther the rotation curve or the gas mass fraction changes 
significantly. 

Based on our analytic solution for the properties of a 
gravitational instability-dominated disk and their depen- 
dence on the gas mass fraction and the infall rate, we are 
able to gain new insight into several processes. We show 
that the velocity dispersions of both the outer H I disks 
of present day galaxies and the main disks of redshift 
~ 2 galaxies can be understood naturally if they are in 
a state of gravitational instability-regulated equilibrium. 
Moreover, we can understand the general progression of 
galactic disks from low values of rotation speed to ve- 
locity dispersion ratio, V max /a, at high redshift to much 
higher values today. This progression is driven primarily 
by a falloff in galaxy accretion rates and secondarily by 
the development of disks with stellar velocity dispersion 
much lager than the gas velocity dispersion, reducing the 
importance of stars in setting the gravitational instabil- 
ity condition. We also predict that the observed range 
of Knax/cr values seen at z ~ 2 is primarily a sequence 
in gas mass fraction. Finally, we use the same model to 
study the velocity dispersions of protostellar disks. We 
show that our results are in good agreement with numer- 
ical simulations of gravitational instability in disks, and 
we predict that velocity dispersions should be transsonic 
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in class and I protostars, dropping to subsonic for class 
II and III sources. 

Although our attention in this paper is focused on cases 
that can be solved analytically or nearly so, we close by 
pointing out that our model, as a result of its grounding 
in the basic equations of fluid dynamics, is also amenable 
to a more general numerical treatment. One can easily 
relax our assumptions of constant gas fraction, negligible 
influence from stellar feedback, and a fixed relationship 
between gas and star velocity dispersion. The result- 
ing equations are identical to the ones we have already 
solved, except that they would need to be solved nu- 
merically. There is no fundamental barrier to doing so 
however, and the result would be a new method for sim- 
ulating the evolution of marginally unstable star-forming 
disks that is intermediate between purely analytic models 
such as those we have pursued here and full numerical 



simulations that can be extremely costly. We plan to 
pursue this avenue in future work. 
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APPENDIX 

A. DERIVATION OF THE TRANSPORT EQUATIONS 

Here we derive the evolution equations for a thin, axisymmetric disk evolving following the general fluid equations 
Jl}, ([2]) and ( |3t including star formation. The d erivation follows the same general outline as the standard treatment 
of disks (e.g. IShul Il992t iBalbus fc Hawlevl 119981 ). with additional terms added to describe star formation and some 
su btleties that aris e in ho w to treat the energy content of supersonic turbulence. We treat these following the method 
of lKrumholz"et"aTI (pOM ) 

Writing out equation (TTJ) in cylindrical coordinates chosen so that the disk lies in the z = plane, dropping terms 
that are zero in axisymmetry, and integrating over z gives equation ([5]), which we repeat here for convenience: 



= --tt^&v) 
at r Or 



t - 1 9 M 
* 2irr dr 



(Al) 



where £ = J pdz is the gas surface density, 
component of the velocity, and we have defined 



J p* dz is the star formation rate per unit area, v r is the radial 

(A2) 



M 



-2nrYiV r 



as the inward radial mass flux. 
Writing out the <fi component of the Navier-Stokes equation ([2]) and performing a similar integration over z yields 



d v r d 



f 1 9 (Jr 
J ^d~r {r 



) dz 



(A3) 



where v$ is the <j) component of the velocity and T rt j, is the rip component of the pressure tensor. Multiplying this 
equation by 2irr 2 gives the evolution equation for the angular momentum ([3]): 

— / 27rr T rct> dz = — 
or / a; 



27rr£ 



di j+Vr dr- J 



-T 



(A4) 



where j — rv^ is the specific angular momentum of the gas and T = J 2irr 2 T r( } > dz is the viscous torque on the gas. 

The gas velocity dispersion is determined by energy conservation. Taking the dot product of v with the Navier-Stokes 
equation @ and using the continuity equation (fT]) to re-arrange, yields 



d (\ 



di\2 PV 



v . v • t - - P y. 



Using the continuity equation (TTJ), we can rewrite the gravitational work term as 

— pv ■ V-ip = —V ■ {pvij;) + ipS7 ■ (pv) 

<9V 



d 

-V ■ (pVTp) - — {pip) + p 



dt 



p*ip- 



Substituting this into equation (|A5[) gives the evolution equation for the non-thermal energy: 



d 



V • pv 



v 



-v-Vp + p- 



v ■ V ■ T — 



v 



(A5) 

(A6) 
(A7) 

(A8) 
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To include the internal energy, we make use ol the hrst law of thermodynamics, equation ([3]). Combining this with 
the continuity equation yields 



^(pe) + V • pv (e+ ~ ) = v • Yi> />,< <I> f V A. 



(A9) 



Adding equations (|A8p and (|A9|) yields the total energy equation: 



( v \ ( iP" v\ dib ( iP 1 

g- t P [ — + e + ipj + V • pv ( y + e + V + M = p^r + v ■ V ■ T + $ - p A — + e + i> ) + T - A. (A10) 



We now integrate over z and use our axisymmetric thin disk assumption to drop all terms involving either z velocities 
or <j) derivatives. This gives 



dt' 



ld_ 

r dr 



= E 



dip 



1 d 



dt 2irr dr 



(QT) 



■4> 



g-c. (ah) 



where (I = v<p/r, Q = jTdz, and C = J Adz. In deriving this equation we have assumed p — £<5(z) and p* = Y;*S(z), 
so in this and all subsequent equations we understand that all the terms in parentheses are to be evaluated in the 
plane z = 0. 

It is convenient to rewrite the kinetic energy plus thermal energy term v 2 /2 + e as a sum of terms representing bulk 
motions on length scales comparable to the radial extent of the galactic disk, small-scale turbulent motions on scales 
comparable to the disk scale height, and true thermal energy: 



1 2 

2 V 



1 / 2 



1 / 2 

5 K 



r 



(A12) 



Here cr n t is the one-dimensional non-thermal velocity dispersion of the turbulent motions on the size scale of the disk 
scale height, of — (2/3)e is the thermal velocity dispersion, and a 1 ~ of + of t EI It is somewhat less clear how to 
evaluate the pressure p in terms of at and <r n t . The microphysical pressure is p = pa\ , but if we are averaging over 
scales much larger than the characteristic size of the turbulent eddies (which is of order the disk scale height), then 
the eddies provide an additional effective pressure, and we will instead have p — pa 2 . Since we are interested in the 
large-scale behavior of disks, we make this microturbulent approximation and adopt p — pa 2 . 

Given this decomposition of the energy, we can simplify equation (|A11|) by dropping small terms. For a rotation- 
dominated disk, v r <C W0. For a disk with dimensionless viscosity a and scale height H, the radial velocity v r ~ 
a(H/r)a. Thus v r <Si a unless the disk is thick {H/r ~ 1) and accretion happens on a dynamical timescale (a ~ 1). 
For this reason we drop the v 2 term. We retain terms of order a com pared to those of order v<p, since we are interested 
in the change in a term of order a. Doing so reduces equation (IA11[) to 



1' 




ld_ 

r dr 



,dip 



1 d 



dt 2irr dr 



(nr)-E, 



-a 2 + ip +0-A (A13) 



where we have rewritten the pressure as p = pa 2 . We can simplify this greatly by using the continuity equation (|A1 1) 
to evaluate the terms on the left-hand side. Doing so, we obtain 



dt 



3a' 



9 < 2 



3a 2 



2i>) 



1 d 



(rT,v r a 2 ) 



1 d 
2ixr dr 



(Q,T) + g-c, 



(A14) 



which is equation (|6|). 



B. SOLUTION OF THE TORQUE EQUATION NEAR THE SINGULARITY AT THE ORIGIN 



Here we obtain the solution to the torque equation (|22|) near x — for constant /3 by means of series expansion. 
Since the nature of the singularity depends on (3, we must handle individual values of /? separately. 



B.l. Flat Rotation Curve ((3 = 0) 
For /3 = the torque equation (|2"2")l reduces to 

\3f g - l) s (3/ s - l)s 2 x 2 \Zf g l) \x) ^ 



5 Note that the coefficient 2/3 in the relation between at and e is 
appropriate for a monatomic ideal gas, i.e. 7 = 5/3; for molecular 
gas 7 can have a different value if the temperature is high enough to 
excite rotational levels of H2 or to induce changes in the ortho- to 



para-H2 ratio. However, for galaxies with a molecule-dominated 
ISM, cr nt is always small compared to a t , so we need not worry 
about a small variations in the coefficient of at ■ 
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We expand r in a power series about x = 0, r = X)^Lo a n% n , to obtain 



(3/ g - I)* 2 



2V2s 3 f gV + oix 
(3/ g - 1)* 2 X 



a 2 [l-2(3/ 5 -l) S 2 ]+5ai/ s ss' 



(3/ fl 



l)s 2 



+ 0{x) = 



Thus the leading coefficients near x — are 



a = 



= -2V2s 3 f g 



0,2 - 



WV2fg 



A J 



1 - 2(3/ s - l)s 2 



(B2) 

(B3) 
(B4) 

(B5) 



For ft 



B.2. Keplerian Rotation Curve (/? = —1/2) 
-1/2 the torque equation ([22]) reduces to 

2 3 / 2 / sS 



(3/ fl -l)s-l(W _, 

(3/g - 1)SX 



4(3/ g - l)s 2 x 3 ' 



3/ fl 

We expand t in a power series about x = 0, r = x 1 ^ 2 X)^o a nX n , to obtain 



1 



X/ a; 



5/2 ' 



-5/2 



3ai 



-3/2 



4(3/ s -l) S 2 x 4(3/ 9 -l) S 2 ' 
Thus the leading coefficients near x = arc 



3ci2 + 2s[5ctos' + 3(3/ g — l)ai] 

4(3/g - 



- 1 / 2 + O(x 1 / 2 ) = 0. 



a = - 
ai=Q 
a2 = 



3Vfg 



(B6) 

(B7) 

(B8) 
(B9) 
(BIO) 

C. NUMERICAL ALGORITHM FOR TIME-DEPENDENT DISKS 

Here we describe our algorithm for numerical solution of the evolution equations (|2"2"|) and (|2"5|) for time-dependent 

disks. Let be the velocity dispersion at time n at the center of cell i, where i runs from 1 to N x . Cell centers are 

located at positions x^ = x^ ^ 1 ^' JVx 1 -' 5 so that the cell spacing in lncc has a uniform value dlax = —(hiXo)/(N x — 1). 

At each time step we obtain the new velocity dispersions using an operator-splitting method in which we treat 

the updates due to the torque explicitly and then perform an implicit diffusion step to suppress spurious numerical 
oscillations. The explicit part of the algorithm, which occurs first in every time step, is: 

1. We compute the spatial derivatives (ds/dx)^ at the center of every grid cell using a minmod slope limiter. 



3 S 
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2. Using s\ n> and (ds/dx)[ , we solve the torque equation (1221) for /3 = and the specified values of x an d f g , 
using the boundary conditions r' = — xq at x = xq and t' = — 1 at x = 1. We solve the equation using the 
method of shooting to a fitting point, with the fitting point chosen in the middle of the computational grid. We 
use an adaptive error control method to maintain accuracy, which is necessary because the torque equation can 
be extremely stiff when \ is small or s is far from the equilibrium solution. This stiffness also limits the smallest 
value of Xq we can use and maintain numerical stability at reasonable computational cost. 



As with s, we evaluate the derivatives 



3. We compute the time derivatives [ds/dty^ in each cell using equation 
of t using a minmod slope limiter. 

4. We set the time step dt = t^ n+1 '> - t^ to dt = 0.02min l [|sf l) /(9s/9x) l ( ™ ) |]. 

5. We set s|"*' > = + dt(ds / dt)[ n \ thereby updating s to time n* for the non-diffusive part of the evolution. 

For the next step, we wish to diffuse the velocity dispersion to prevent the development of grid-scale numerical 
oscillations. In order to guarantee energy conservation, we diffuse the kinetic energy rather than diffusing s directly. 
We define the dimensionless kinetic energy per unit area in a computational cell by 



14 



and we evolve this following 



dk 



KdiffV k 



K diff d 2 k 
x 2 d\nx' 



(C2) 



where Kas is the diffusion coefficient, and in the second step we have evaluated the V 2 operator on our cylindrical 
logarithmic grid. Provided that we set dk/dhix = at our inner and outer boundaries, evolution under this equation 
does not change the total amount of kinetic energy on the grid. We discretize equation (|C2p using centered spatial 
differences and fully implicit temporal differences: 



(n.) 



Kdiff 



ill 



>+!) _l 



2k 



(n+l) 



(din a;) 2 



(C3) 



with the boundary conditions that /cq™ +1 ^ = k^ +1 ^ and k^ + }\ = kpf + ■ Rewriting this in matrix form, we have 



where k(" +1 ) and 



M . k ("+1) = k (n,) ; ( C4 ) 

are the vectors of ki values at times n+l and n» , respectively, and the matrix M has elements 



Mi. 



Kdifidt 
(din a:) 2 



X 



(C5) 



N. r 



Since M is tridiagonal matrix, it is easy to solve equation (IC4[) exactly. Thus to take our diffusion step we simply 
compute k^™*) from s^ n *' (equation lClj) . solve equation ()C4[) for k^ n+1 -', and use this to compute s(™ +1 - 1 , thus completing 
the timestep. We find that Kdiff = 0.005 (for \ — 0.01) or Kdiff = 0.01 (for x = 0.1) is sufficient to suppress numerical 
oscillations without significantly changing the solution. 
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